* dependencies:
*   estout
*   reghdfejl and julia
*   blindschemes
*   grc1leg2


cap cd "D:\OneDrive - Open Philanthropy Project\Vaccines\Summan, Nandi, and Bloom"
cap cd "/Users/davidroodman/Library/CloudStorage/OneDrive-OpenPhilanthropyProject/Vaccines/Summan, Nandi, and Bloom"

graph set window fontface "LM Roman 9"
set scheme plotplain // plottig

***
*** Prep data
***
{
//   cap noi odbc load, clear dsn("SNB") exec("select * from [Analysis data]")
//   if !_rc {
//     compress
//     saveold "Analysis data", version(12) replace
//   }
  use "Analysis data", clear

  recode General_Education* headedu* (1/7 = 1)  // http://microdata.gov.in/nada43/index.php/catalog/127/datafile/F4/V161
  recode Relation_to_Head* (5 = 3) (. 4 8 9 = 0)  // http://microdata.gov.in/nada43/index.php/catalog/127/datafile/F4/V157
  gen byte married = 2.Marital_Status

  bysort HHID (PID): gen headedu    = General_Education[1]
  by     HHID      : gen headage    = Age[1]
  by     HHID      : gen femalehead = !Male[1]
  recode headedu (. = 0)

  preserve
  use "2023-05-18 SNB reply/nfhs4_migrant_clean", clear
  probit migrant Age Agesq nevermarried child grandchild mpceq2 mpceq3 mpceq4 mpceq5 female head wifehead headage femalehead rural sc st obc muslim christian sikh HH_Size if work & birthyear>=1985 & birthyear<=1990
  restore
  gen double migrant85_90SNB = normprob(            _b[female]*femaleSNB +_b[Age]*AgeSNB +_b[Agesq]*AgeSNB^2	+_b[nevermarried]*nevermarriedSNB	+_b[head]*headSNB	+_b[wifehead]*wifeheadSNB	+_b[child]*childSNB	+_b[grandchild]*grandchildSNB	+_b[headage]*headageSNB	+_b[femalehead]*femaleheadSNB	+_b[rural]*ruralSNB	+_b[sc]*scSNB	+_b[st]*stSNB	+_b[obc]*obcSNB	+_b[muslim]*muslimSNB	+_b[christian]*christianSNB	+_b[sikh]*sikhSNB	+_b[HH_Size]*HH_SizeSNB	+_b[mpceq2]*mpceq2SNB	+_b[mpceq3]*mpceq3SNB	+_b[mpceq4]*mpceq4SNB	+_b[mpceq5]*mpceq5SNB)
  gen double migrant85_90    = normprob(_b[_cons] + _b[female]*femaleSNB +_b[Age]*AgeSNB +_b[Agesq]*AgeSNB^2	+_b[nevermarried]*nevermarriedSNB	+_b[head]*headSNB	+_b[wifehead]*wifeheadSNB	+_b[child]*childSNB	+_b[grandchild]*grandchildSNB	+_b[headage]*headageSNB	+_b[femalehead]*femaleheadSNB	+_b[rural]*ruralSNB	+_b[sc]*scSNB	+_b[st]*stSNB	+_b[obc]*obcSNB	+_b[muslim]*muslimSNB	+_b[christian]*christianSNB	+_b[sikh]*sikhSNB	+_b[HH_Size]*HH_SizeSNB	+_b[mpceq2]*mpceq2SNB	+_b[mpceq3]*mpceq3SNB	+_b[mpceq4]*mpceq4SNB	+_b[mpceq5]*mpceq5SNB)

  replace Date_of_Survey = date(string(Date_of_Survey, "%06.0f"), "DM20Y")
  replace Date_of_Survey = Date_of_Survey - 366 if Date_of_Survey > td(30jun2012)  // Amit Summan suggestion: pull a few post-survey observations back a year
  format %td Date_of_Survey
  gen int Month_of_Survey = mofd(Date_of_Survey)

  gen svyyear = yofd(Date_of_Survey)
  gen int birthyear    = svyyear    - Age
  gen int birthyearSNB = svyyearSNB - AgeSNB
  gen byte TSNB = birthyearSNB >= uipSNB
  replace uip = round(uip, 1)
  recast int uip
  replace uip = . if uip_SD > .1
  gen byte T = birthyear >= uip

  gen double lwage = ln(Wages) if !inlist(Usual_Principal_Activity_Status, 11,12,21)  // don't count self-employment wages
  gen double ldwage = lwage - ln(Days_worked)  // convert to wages per time worked

  * compute CPIs as function of dates, interpolating from monthly data, 6/2011-12/2012, from https://data.imf.org/?sk=fffaba71-e02b-4e2c-87bc-45cf050eb5bc&hide_uv=1
  gen CPI = .
  gen byte hasdate = svyyear < .
  mata {
    CPI = 107.44 \ 109.71 \ 110.28 \ 111.98 \ 112.55 \ 113.12 \ 112.63 \ 113.02 \ 113.62 \ 114.51 \ 116.10 \ 117.19 \ 118.57 \ 120.36 \ 121.85 \ 122.94 \ 123.83 \ 124.32 \ 124.52
    dCPI = CPI[|2\.|] - CPI[|.\length(CPI)-1|]
    st_view(surveyCPI=., ., "CPI", "hasdate")
    date = st_data(., "Date_of_Survey", "hasdate")
    y = year(date); m = month(date); d = day(date); my = mofd(date)
    i = y * 12 + m :- (2011*12+6)  // indexes in CPI table for last value before date
    surveyCPI[,] = CPI[i] + d:/(dofm(my:+1) - dofm(my)) :* dCPI[i]  // interpolate
  }

  gen double lwagereal = lwage - ln(CPI/100)
  gen double lmpce = ln(mpce)
  gen double lmpcereal = lmpce - ln(CPI/100)
  
  sum WEIGHT, detail
  gen WTtr = min(WEIGHT, r(p50)+5*(r(p75)-r(p50)))  // clip extreme weights to median + 5 * IQR (Potter and Zheng 2015)
}


global controlsSNB urban Male marriedSNB HH_Size headage i.Social_Group i(2/4)bn.Religion i(1 2 3 6)bn.Relation_to_Head femalehead i(8 10 12 13)bn.(General_Education headedu) migrant85_90SNB
global controls    urban Male married    HH_Size headage i.Social_Group i(2/4)bn.Religion i(1 2 3 6)bn.Relation_to_Head femalehead i(8 10 12 13)bn.(General_Education headedu) migrant85_90

* best replication; then monthly survey effects instead; then inflation adjusted
cap erase esttab.rtf
local keepsing keepsing
foreach SNB in SNB "" {
  foreach depvar in lwage lmpce {
    eststo clear
    eststo: reghdfejl `depvar'     T`SNB' $controls`SNB' if inrange(birthyear`SNB',1985,1990)          , cluster(District_Code) a(Age#District_Code) `keepsing'  // matches original sample size & analysis.do, line 128, sent 5/18/23 
    eststo: reghdfejl `depvar'     T`SNB' $controls`SNB' if inrange(birthyear`SNB',1985,1990) [aw=WTtr], cluster(District_Code) a(Age#District_Code) `keepsing'  // matches original sample size & analysis.do, line 128, sent 5/18/23 
    eststo: reghdfejl `depvar'     T`SNB' $controls`SNB' if inrange(birthyear`SNB',1985,1990)          , cluster(District_Code) a(Age#District_Code Month_of_Survey) `keepsing'  // matches original sample size & analysis.do, line 128, sent 5/18/23 
    eststo: reghdfejl `depvar'     T`SNB' $controls`SNB' if inrange(birthyear`SNB',1985,1990) [aw=WTtr], cluster(District_Code) a(Age#District_Code Month_of_Survey) `keepsing'  // matches original sample size & analysis.do, line 128, sent 5/18/23 
    eststo: reghdfejl `depvar'real T`SNB' $controls`SNB' if inrange(birthyear`SNB',1985,1990)          , cluster(District_Code) a(Age#District_Code) `keepsing'  // matches original sample size & analysis.do, line 128, sent 5/18/23 
    eststo: reghdfejl `depvar'real T`SNB' $controls`SNB' if inrange(birthyear`SNB',1985,1990) [aw=WTtr], cluster(District_Code) a(Age#District_Code) `keepsing'  // matches original sample size & analysis.do, line 128, sent 5/18/23 
    esttab using esttab.rtf, append b(3) se(a2) nocons mtitles rename(T`SNB' T) keep(T) stat(N, label("Observations") fmt(%7.0fc)) nonumber nonote nogap label nostar noline msign("–") fonttbl(\f0\fnil LM Roman 9;)
  }
  local keepsing
}

* wages and consumption rise during the survey
cap drop at
gen at = td(1jul2011) in 1
forvalues d=1/4 {
  local depvar: word `d' of TSNB lwage lmpce urban
  local depname: word `d' of Treated "Log wages" "Log expenditure per person" Urban
  local ylabs  : word `d' of "" 7(.2)7.4 7.2(.05)7.35
  local xscale: word `d' of off
  cap drop fit1
  lpoly `depvar' Date_of_Survey if inrange(birthyearSNB,1985,1990), bw(30) gen(fit1) at(at) nograph
  cap drop fit2
  lpoly `depvar' Date_of_Survey if inrange(birthyearSNB,1975,1980), bw(30) gen(fit2) at(at) nograph

  twoway lpoly `depvar' Date_of_Survey if inrange(birthyearSNB,1985,1990), bw(30) lstyle(p1) || ///
         lfit  `depvar' Date_of_Survey if inrange(birthyearSNB,1985,1990), lstyle(p3) || ///
         lpoly `depvar' Date_of_Survey if inrange(birthyearSNB,1975,1980), bw(30) lstyle(p1) || ///
         lfit  `depvar' Date_of_Survey if inrange(birthyearSNB,1975,1980), lstyle(p3) text(`=fit2[1]' `=at[1]' "Born" "1975–80", place(ne) size(small) just(left) margin(0 0 1 0)) ///
         text(`=fit1[1]' `=at[1]' "Born" "1985–90", place(ne) size(small) just(left) margin(0 0 1 0)) ///
         title("`depname'") ylab(`ylabs',format(%3.2f)) ///
         xlab(`=td(1jul2011)'(183)`=td(1jul2012)', format(%tdMon-ccyy)) legend(off) ///
         name(`depvar'trend, replace) 
}
graph combine TSNBtrend lwagetrend lmpcetrend, rows(1) iscale(1) imargin(0 3 3 0) b1title(Survey date, margin(medsmall)) 
graph export trends.png, replace width(3000)
graph export trends.tif, replace width(3000)

* Checks for slope differences between young and old
cap drop old
gen byte old = birthyearSNB < 1985
gen _Date_of_Survey = Date_of_Survey / 366
reg lwage i.old old#c._Date_of_Survey if inrange(birthyearSNB,1985,1990) | inrange(birthyearSNB,1975,1980), cluster(District_Code)
reg lmpce i.old old#c._Date_of_Survey if inrange(birthyearSNB,1985,1990) | inrange(birthyearSNB,1975,1980), cluster(District_Code)


***
*** Randomization inference: scramble the treatment variable; two-tail test of null that point estimate = mean of empirical distribution
***
{
cap program drop putuip
program define putuip
  preserve
  collapse (firstnm) `1', by(clustid)
  jl PutVarsToMat `1', dest(uip)
  jl: uip = Int16.(uip);
  restore
end


jl: using Distributed; rmprocs(workers()); addprocs(6)
jl: @everywhere using FixedEffectModels, Random, SharedArrays

global M 100
jl: b = SharedVector{Float64}($M);

recode uip* (. = 2100)  // treat missing treatment year as treatment coming much later; recode so Julia doesn't have to deal with missings

foreach SNB in SNB /*""*/ {
  local controls HH_Size headage i(2/4)bn.Religion i(1 2 3 6)bn.Relation_to_Head i(8 10 12 13)bn.(General_Education headedu) migrant85_90`SNB'
  local FE urban Male married`SNB' Social_Group femalehead
  local FEexp fe(urban)+fe(Male)+fe(married`SNB')+fe(Social_Group)+fe(femalehead)

  forvalues w=1/2 {
    local wtexp  : word `w' of "" [aw=WTtr]
    local wtexpjl: wor* dependencies:
*   estout
*   reghdfejl and julia
*   blindschemes
*   grc1leg2


cap cd "D:\OneDrive - Open Philanthropy Project\Vaccines\Summan, Nandi, and Bloom"
cap cd "/Users/davidroodman/Library/CloudStorage/OneDrive-OpenPhilanthropyProject/Vaccines/Summan, Nandi, and Bloom"

graph set window fontface "LM Roman 9"
set scheme plotplain // plottig

***
*** Prep data
***
{
//   cap noi odbc load, clear dsn("SNB") exec("select * from [Analysis data]")
//   if !_rc {
//     compress
//     saveold "Analysis data", version(12) replace
//   }
  use "Analysis data", clear

  recode General_Education* headedu* (1/7 = 1)  // http://microdata.gov.in/nada43/index.php/catalog/127/datafile/F4/V161
  recode Relation_to_Head* (5 = 3) (. 4 8 9 = 0)  // http://microdata.gov.in/nada43/index.php/catalog/127/datafile/F4/V157
  gen byte married = 2.Marital_Status

  bysort HHID (PID): gen headedu    = General_Education[1]
  by     HHID      : gen headage    = Age[1]
  by     HHID      : gen femalehead = !Male[1]
  recode headedu (. = 0)

  preserve
  use "2023-05-18 SNB reply/nfhs4_migrant_clean", clear
  probit migrant Age Agesq nevermarried child grandchild mpceq2 mpceq3 mpceq4 mpceq5 female head wifehead headage femalehead rural sc st obc muslim christian sikh HH_Size if work & birthyear>=1985 & birthyear<=1990
  restore
  gen double migrant85_90SNB = normprob(            _b[female]*femaleSNB +_b[Age]*AgeSNB +_b[Agesq]*AgeSNB^2	+_b[nevermarried]*nevermarriedSNB	+_b[head]*headSNB	+_b[wifehead]*wifeheadSNB	+_b[child]*childSNB	+_b[grandchild]*grandchildSNB	+_b[headage]*headageSNB	+_b[femalehead]*femaleheadSNB	+_b[rural]*ruralSNB	+_b[sc]*scSNB	+_b[st]*stSNB	+_b[obc]*obcSNB	+_b[muslim]*muslimSNB	+_b[christian]*christianSNB	+_b[sikh]*sikhSNB	+_b[HH_Size]*HH_SizeSNB	+_b[mpceq2]*mpceq2SNB	+_b[mpceq3]*mpceq3SNB	+_b[mpceq4]*mpceq4SNB	+_b[mpceq5]*mpceq5SNB)
  gen double migrant85_90    = normprob(_b[_cons] + _b[female]*femaleSNB +_b[Age]*AgeSNB +_b[Agesq]*AgeSNB^2	+_b[nevermarried]*nevermarriedSNB	+_b[head]*headSNB	+_b[wifehead]*wifeheadSNB	+_b[child]*childSNB	+_b[grandchild]*grandchildSNB	+_b[headage]*headageSNB	+_b[femalehead]*femaleheadSNB	+_b[rural]*ruralSNB	+_b[sc]*scSNB	+_b[st]*stSNB	+_b[obc]*obcSNB	+_b[muslim]*muslimSNB	+_b[christian]*christianSNB	+_b[sikh]*sikhSNB	+_b[HH_Size]*HH_SizeSNB	+_b[mpceq2]*mpceq2SNB	+_b[mpceq3]*mpceq3SNB	+_b[mpceq4]*mpceq4SNB	+_b[mpceq5]*mpceq5SNB)

  replace Date_of_Survey = date(string(Date_of_Survey, "%06.0f"), "DM20Y")
  replace Date_of_Survey = Date_of_Survey - 366 if Date_of_Survey > td(30jun2012)  // Amit Summan suggestion: pull a few post-survey observations back a year
  format %td Date_of_Survey
  gen int Month_of_Survey = mofd(Date_of_Survey)

  gen svyyear = yofd(Date_of_Survey)
  gen int birthyear    = svyyear    - Age
  gen int birthyearSNB = svyyearSNB - AgeSNB
  gen byte TSNB = birthyearSNB >= uipSNB
  replace uip = round(uip, 1)
  recast int uip
  replace uip = . if uip_SD > .1
  gen byte T = birthyear >= uip

  gen double lwage = ln(Wages) if !inlist(Usual_Principal_Activity_Status, 11,12,21)  // don't count self-employment wages
  gen double ldwage = lwage - ln(Days_worked)  // convert to wages per time worked

  * compute CPIs as function of dates, interpolating from monthly data, 6/2011-12/2012, from https://data.imf.org/?sk=fffaba71-e02b-4e2c-87bc-45cf050eb5bc&hide_uv=1
  gen CPI = .
  gen byte hasdate = svyyear < .
  mata {
    CPI = 107.44 \ 109.71 \ 110.28 \ 111.98 \ 112.55 \ 113.12 \ 112.63 \ 113.02 \ 113.62 \ 114.51 \ 116.10 \ 117.19 \ 118.57 \ 120.36 \ 121.85 \ 122.94 \ 123.83 \ 124.32 \ 124.52
    dCPI = CPI[|2\.|] - CPI[|.\length(CPI)-1|]
    st_view(surveyCPI=., ., "CPI", "hasdate")
    date = st_data(., "Date_of_Survey", "hasdate")
    y = year(date); m = month(date); d = day(date); my = mofd(date)
    i = y * 12 + m :- (2011*12+6)  // indexes in CPI table for last value before date
    surveyCPI[,] = CPI[i] + d:/(dofm(my:+1) - dofm(my)) :* dCPI[i]  // interpolate
  }

  gen double lwagereal = lwage - ln(CPI/100)
  gen double lmpce = ln(mpce)
  gen double lmpcereal = lmpce - ln(CPI/100)
  
  sum WEIGHT, detail
  gen WTtr = min(WEIGHT, r(p50)+5*(r(p75)-r(p50)))  // clip extreme weights to median + 5 * IQR (Potter and Zheng 2015)
}


global controlsSNB urban Male marriedSNB HH_Size headage i.Social_Group i(2/4)bn.Religion i(1 2 3 6)bn.Relation_to_Head femalehead i(8 10 12 13)bn.(General_Education headedu) migrant85_90SNB
global controls    urban Male married    HH_Size headage i.Social_Group i(2/4)bn.Religion i(1 2 3 6)bn.Relation_to_Head femalehead i(8 10 12 13)bn.(General_Education headedu) migrant85_90

* best replication; then monthly survey effects instead; then inflation adjusted
cap erase esttab.rtf
local keepsing keepsing
foreach SNB in SNB "" {
  foreach depvar in lwage lmpce {
    eststo clear
    eststo: reghdfejl `depvar'     T`SNB' $controls`SNB' if inrange(birthyear`SNB',1985,1990)          , cluster(District_Code) a(Age#District_Code) `keepsing'  // matches original sample size & analysis.do, line 128, sent 5/18/23 
    eststo: reghdfejl `depvar'     T`SNB' $controls`SNB' if inrange(birthyear`SNB',1985,1990) [aw=WTtr], cluster(District_Code) a(Age#District_Code) `keepsing'  // matches original sample size & analysis.do, line 128, sent 5/18/23 
    eststo: reghdfejl `depvar'     T`SNB' $controls`SNB' if inrange(birthyear`SNB',1985,1990)          , cluster(District_Code) a(Age#District_Code Month_of_Survey) `keepsing'  // matches original sample size & analysis.do, line 128, sent 5/18/23 
    eststo: reghdfejl `depvar'     T`SNB' $controls`SNB' if inrange(birthyear`SNB',1985,1990) [aw=WTtr], cluster(District_Code) a(Age#District_Code Month_of_Survey) `keepsing'  // matches original sample size & analysis.do, line 128, sent 5/18/23 
    eststo: reghdfejl `depvar'real T`SNB' $controls`SNB' if inrange(birthyear`SNB',1985,1990)          , cluster(District_Code) a(Age#District_Code) `keepsing'  // matches original sample size & analysis.do, line 128, sent 5/18/23 
    eststo: reghdfejl `depvar'real T`SNB' $controls`SNB' if inrange(birthyear`SNB',1985,1990) [aw=WTtr], cluster(District_Code) a(Age#District_Code) `keepsing'  // matches original sample size & analysis.do, line 128, sent 5/18/23 
    esttab using esttab.rtf, append b(3) se(a2) nocons mtitles rename(T`SNB' T) keep(T) stat(N, label("Observations") fmt(%7.0fc)) nonumber nonote nogap label nostar noline msign("–") fonttbl(\f0\fnil LM Roman 9;)
  }
  local keepsing
}

* wages and consumption rise during the survey
cap drop at
gen at = td(1jul2011) in 1
forvalues d=1/4 {
  local depvar: word `d' of TSNB lwage lmpce urban
  local depname: word `d' of Treated "Log wages" "Log expenditure per person" Urban
  local ylabs  : word `d' of "" 7(.2)7.4 7.2(.05)7.35
  local xscale: word `d' of off
  cap drop fit1
  lpoly `depvar' Date_of_Survey if inrange(birthyearSNB,1985,1990), bw(30) gen(fit1) at(at) nograph
  cap drop fit2
  lpoly `depvar' Date_of_Survey if inrange(birthyearSNB,1975,1980), bw(30) gen(fit2) at(at) nograph

  twoway lpoly `depvar' Date_of_Survey if inrange(birthyearSNB,1985,1990), bw(30) lstyle(p1) || ///
         lfit  `depvar' Date_of_Survey if inrange(birthyearSNB,1985,1990), lstyle(p3) || ///
         lpoly `depvar' Date_of_Survey if inrange(birthyearSNB,1975,1980), bw(30) lstyle(p1) || ///
         lfit  `depvar' Date_of_Survey if inrange(birthyearSNB,1975,1980), lstyle(p3) text(`=fit2[1]' `=at[1]' "Born" "1975–80", place(ne) size(small) just(left) margin(0 0 1 0)) ///
         text(`=fit1[1]' `=at[1]' "Born" "1985–90", place(ne) size(small) just(left) margin(0 0 1 0)) ///
         title("`depname'") ylab(`ylabs',format(%3.2f)) ///
         xlab(`=td(1jul2011)'(183)`=td(1jul2012)', format(%tdMon-ccyy)) legend(off) ///
         name(`depvar'trend, replace) 
}
graph combine TSNBtrend lwagetrend lmpcetrend, rows(1) iscale(1) imargin(0 3 3 0) b1title(Survey date, margin(medsmall)) 
graph export trends.png, replace width(3000)
graph export trends.tif, replace width(3000)

* Checks for slope differences between young and old
cap drop old
gen byte old = birthyearSNB < 1985
gen _Date_of_Survey = Date_of_Survey / 366
reg lwage i.old old#c._Date_of_Survey if inrange(birthyearSNB,1985,1990) | inrange(birthyearSNB,1975,1980), cluster(District_Code)
reg lmpce i.old old#c._Date_of_Survey if inrange(birthyearSNB,1985,1990) | inrange(birthyearSNB,1975,1980), cluster(District_Code)


***
*** Randomization inference: scramble the treatment variable; two-tail test of null that point estimate = mean of empirical distribution
***
{
cap program drop putuip
program define putuip
  preserve
  collapse (firstnm) `1', by(clustid)
  jl PutVarsToMat `1', dest(uip)
  jl: uip = Int16.(uip);
  restore
end


jl: using Distributed; rmprocs(workers()); addprocs(6)
jl: @everywhere using FixedEffectModels, Random, SharedArrays

global M 100000
jl: b = SharedVector{Float64}($M);

recode uip* (. = 2100)  // treat missing treatment year as treatment coming much later; recode so Julia doesn't have to deal with missings

foreach SNB in SNB "" {
  local controls HH_Size headage i(2/4)bn.Religion i(1 2 3 6)bn.Relation_to_Head i(8 10 12 13)bn.(General_Education headedu) migrant85_90`SNB'
  local FE urban Male married`SNB' Social_Group femalehead
  local FEexp fe(urban)+fe(Male)+fe(married`SNB')+fe(Social_Group)+fe(femalehead)

  forvalues w=1/2 {
    local wtexp  : word `w' of "" [aw=WTtr]
    local wtexpjl: word `w' of "" ", weights=:WTtr"

    forvalues d=1/2 {
      local depvar : word `d' of lwage lmpce
      local depname: word `d' of "Log wage" "Log per-capita expenditure"

      reghdfejl `depvar' T`SNB' `controls' `wtexp' if inrange(birthyear`SNB',1985,1990), a(`FE' Age`SNB'#District_Code)
      local b0 = _b[T`SNB']

      preserve

      keep if e(sample)
      egen int clustid = group(District_Code)
      fvrevar `controls'
      local _controls `r(varlist)'
      local controlsjl: subinstr local _controls " " "+", all

      putuip uip`SNB'    
      jl PutVarsToDF `depvar' T`SNB' Age`SNB' birthyear`SNB' clustid `_controls' `FE' WTtr, doubleonly nomissing
      jl: df.birthyear`SNB' = Vector{Int16}(df.birthyear`SNB');
      jl: df.clustid = Vector{Int16}(df.clustid);
      
      // For speed partial out FE from non-randomized variables just once
      jl: dfp = partial_out(df, @formula(`depvar' + `controlsjl' ~ `FEexp' + fe(Age`SNB')&fe(clustid)) `wtexpjl')[1];
      jl: df[!,names(dfp)] = dfp;  // overwrite original vars with partialled-out versions

      jl: partialf = @formula(T`SNB' ~ `FEexp' + fe(Age`SNB')&fe(clustid));  // regression specification for partialling FE out of randomized dep var
      jl: f = @formula(`depvar' ~ 0 + T`SNB' + `controlsjl');  // specification for primary regression

      jl: @sync @distributed for m in 1:$M                                           ///
            df.T`SNB' .= df.birthyear`SNB' .>= getindex.(Ref(shuffle(uip)), df.clustid); ///
            df.T`SNB' = partial_out(df, partialf `wtexpjl')[1].T`SNB';                   ///
            b[m] = reg(df, f `wtexpjl', nthreads=1, drop_singletons=false, progress_bar=false).coef[1]                          ///
          end

      drop _all
      set obs $M
      jl GetVarsFromMat b, source(b)
      sum b, meanonly
      local mean = r(mean)
      qui count if b<. & abs(`b0' - r(mean)) < abs(b - `mean')  // two-tailed mass
      twoway kdensity b, lcolor(gs10) xline(`b0', lcolor(gs2) lpat(solid)) xline(`mean', lcolor(gs10) lpat(dash)) note("") name(`depvar'`w', replace) xtitle("") yscale(off) ///
        title(`"`depname', `:word `w' of unweighted weighted' regression"', size(small)) ///
        ylab(,nogrid) ///
        xlab(0, custom `=cond(`w'==2,"","nolab notick grid")') xlab(`b0' `" "`:display %4.3f `b0''" "({it:p} = `:display %3.2f r(N)/$M')" "', add custom labgap(10pt) labcolor(gs2) tlcolor(gs2)) ///
        xlab(`mean' `" "`:display %4.3f `mean''" "', add custom labcolor(gs10) tlcolor(gs10)) xscale(extend) || ///
        line b b b in 1, lcolor(gs10 gs2) lpat(dash solid) legend(order(2 3) lab(2 "Randomization mean") lab(3 "Point estimate"))  // fake plot to make legend
        
      restore
    }
  }
  grc1leg2 lwage1 lmpce1 lwage2 lmpce2, iscale(1) imargin(1 1 3 3) lrows(1) legscale(medsmall) xcommon
  graph export RI`SNB'.png, replace width(3000)
  graph export RI`SNB'.tif, replace width(3000)
  }
}